#ifndef DENSITY_HPP
#define DENSITY_HPP

#include <armadillo>
#include <boost/math/distributions/normal.hpp>

using namespace arma;


namespace Density{

// Abstract class for tempering (see implementation of Likelihood lower)

class Tempering{
	public:
		Tempering(double *Y,mat X){
			_Y=Y;
			_X=X;
			_p=X.n_cols;
			_n=X.n_rows;
			_b=0;
			_phi_n=0;
		}
		virtual double Likelihood(mat)=0;
		virtual double Prior(mat)=0;
		
		inline double Density_utn(mat theta, int n){
			return (this->Likelihood(theta))*(_phi_n);
		}
		inline void Weight(mat theta,double *W,int n)
		{
			double temp=0;
			double bar=_phi_n-_phi_n1;
			//cout << "2 " << bar << " 3";
			int M=theta.n_cols;
			double *foo=new double[M];
			for(int i=0;i<M;i++)
			{	temp=(this->Likelihood(theta.col(i)));
				foo[i]=temp*bar;	
				W[i]*=exp(foo[i]);
			//	cout << W[i] << "\n";
			}

			delete[] foo;
		}
		inline double* Eta_star(double phi, double *W, mat theta)
		{
			double temp=0;
			double bar=phi;//-_phi_n1;
			//cout << "2 " << bar << " 3";
			int M=theta.n_cols;
			double *foo=new double[M];// A MODIF
			for(int i=0;i<M;i++)
			{	temp=(this->Likelihood(theta.col(i).t()));
				foo[i]=W[i]*exp(temp*bar);	
			}
			return foo;
		}
		inline void Set_Phi(double phin){
			_phi_n1=_phi_n;
			_phi_n=phin;
			
		}

		inline void Set_p(int p){ _p=p;}
		inline int Get_p(void){return _p;}
		inline double Get_bn(void){return _phi_n;}
		inline int Get_n(void){return _n;}		
		inline	mat Get_X(void){return _X;}
		inline	double* Get_Y(void){return _Y;}
	private:
		double _phi_n;
		double _phi_n1;
		mat _X;
		double *_Y;
		int _p;
		int _b;
		int _n;
};

////////////////////////////////////////////////////////////////////////////////
// Probit IBIS, Should have code of "General" IBIS here
class IBIS{
	public:
		IBIS(double *Y,mat X){
			_Y=Y;
			_X=X;
			_p=X.n_cols;
			_n=X.n_rows;
			_b=0;
		}
		virtual double	Density(mat theta, int n, int p)=0;//Density from n to n+p==> virtual pure, (class becomes Density_IBIS)+ add implementation elsewhere of the specific function....

		virtual double Prior(mat)=0;
		virtual mat Get_v(void)=0;

	//from begining
		inline double 	Density_utn( mat theta, int n)
		{
			return	Density(theta,0,n)+(Prior(theta));
		}
			
		inline void  Weight(mat theta, double *W, int n)
		{
			int M=theta.n_cols;
			for(int i=0;i<M;i++)
			{
				W[i]*=exp(Density(theta.col(i),n,_b));
			}

		}
		inline	int Get_p(void){return _p;}
		inline	int 	Get_n(void){return _n;}
		inline	mat Get_X(void){return _X;}
		inline	double* Get_Y(void){return _Y;}
		inline  double Get_bn(void){return 1;}
		inline  void Set_X(int i){
			mat X(_n,_p);
			double *Y= new double[_n]; 
			int k=0;
			for( int j=0;j<_n;j++)
			{
				if(j!=i)
				{
					Y[k]=_Y[i];
					X(k,span::all)=_X(k,span::all);
					k++;
				}else{
					Y[_n-1]=_Y[i];
					X(_n-1,span::all)=_X(_n-1,span::all);
				}				
			}
			for( int j=0;j<_n;j++)
			{
				_Y[j]=Y[j];
			}
			delete[] Y;
			_X=X;

		}
	private:
			mat _X;
			double *_Y;
			int _p;
			int _b;
			int _n;
};


// Abstract class for tempering (see implementation of Likelihood lower)

class GeomBridge{
	public:
		GeomBridge(double *Y,mat X,Distribution::Distribution *S){
			_Y=Y;
			_X=X;
			_p=X.n_cols;
			_n=X.n_rows;
			_b=0;
			_phi_n=0;
			_Start=S;
		}
		virtual double Likelihood(mat)=0;
		virtual double Prior(mat)=0;
		virtual mat GradLik(mat)=0;
		
		inline double Density_utn(mat theta, int n){
			return ((this->Likelihood(theta))+this->Prior(theta))*(_phi_n)+(1-_phi_n)*(_Start->d(theta,1,0));
		}
		mat Grad(mat theta){
			return (GradLik(theta))*(_phi_n)-(1-_phi_n)*inv(_Start->Get_s())*(theta-_Start->Get_m());
		}
		inline double Weight(mat theta,double *W,int n, int Bool)
		{
			double temp=0;
			double bar=_phi_n-_phi_n1;
			double temp2=0;
			double temp3=0;
			//cout << "2 " << bar << " 3";
			int M=theta.n_cols;
			if(Bool==1)
			{
				
				double *foo=new double[M];
				for(int i=0;i<M;i++)
				{	
					temp=(this->Likelihood(theta.col(i)))+this->Prior(theta.col(i))-_Start->d(theta.col(i),1,0);
					W[i]+=(temp)*bar;	
					//cout << W[i] << "\n";
				}
				delete[] foo;
			}else{
				for(int i=0;i<M;i++)
				{	
					temp=_Targ(i,0);
					W[i]+=(temp)*bar;	
					//cout << W[i] << "\n";
				}
			}
		}
		inline double* Eta_star(double phi, double *W, mat theta,int Bool)
		{
			double temp=0;
			double temp0=0;
			double temp2=0;
			double temp3=0;
			double bar=phi;//-_phi_n1;
			//cout << "2" << bar << "3";
			int M=theta.n_cols;
			double *foo=new double[M];
			if(Bool==1)
			{
				mat Temp(M,1);
				for(int i=0;i<M;i++)
				{	
				//	cout << theta.col(i);
					Temp(i,0)=(this->Likelihood(theta.col(i)))+(this->Prior(theta.col(i)))-_Start->d(theta.col(i),1,0);
					foo[i]=W[i]+(Temp(i,0))*bar;	
					
				/*	
					temp0=(this->Likelihood(theta.col(0)));
					temp2=(this->Prior(theta.col(0)));
					temp3=-_Start->d(theta.col(0),1,0);
					cout << "temp2 " << temp2 <<"\n";
					cout << "temp0 " << temp0 <<"\n";
					cout << "temp3 " << temp3 << " "<< as_scalar(sum(sum(theta.col(i))))<<"\n";
				//	cout << "foo: "<< foo[i]  << " ";
					foo[i]=W[i]+bar*(temp0+temp2+temp3);//(Temp(i,0))*bar;	
				*/	
				}
				_Targ=Temp;
			}else{
				for(int i=0;i<M;i++)
				{	
					foo[i]=W[i]+(_Targ(i,0))*bar;	
					//cout << "foo: "<< temp  << " ";
				}

			}
			//cout << foo <<"\n";
			return foo;
		}
		inline void Set_Phi(double phin){
			_phi_n1=_phi_n;
			_phi_n=phin;
			
		}
		inline int Set_p(int p){ _p=p;}
		inline int Get_p(void){return _p;}
		inline double Get_bn(void){return _phi_n;}
		inline int Get_n(void){return _n;}		
		inline	mat Get_X(void){return _X;}
		inline	double Get_phi_n(void){return _phi_n;}
		inline	double* Get_Y(void){return _Y;}
		inline Distribution::Distribution* Get_Start(){return _Start;}
	private:
		double _phi_n;
		double _phi_n1;
		mat _X;
		double *_Y;
		int _p;
		int _b;
		int _n;
		mat _Targ;
		Distribution::Distribution *_Start;
};

}//End namespace
#endif

